Carbon Ignition in Type la Supernovae: An Analytic Model 



m 
o 
o 

(N 

< 



S. E. Woosley 

Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064 

woosleyOucolick . org 

S. Wunsch 

Combustion Research Facility 
Sandia National Laboratories, Livermore, CA 94551-0969 
sewunscSca. SEindia.gov 
and 



(N 
> 

in 
\o 
in 
i> 
o 
cn 
o 

:^ 

o 
U 



X 



M. Kuhlen 

Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064 

mqkOucolick . org 



ABSTRACT 

The observable properties of a Type la supernova are sensitive to how the nuclear runaway 
ignites in a Chandrasekhar mass white dwarf - at a single point at its center, off-center, or at 
multiple points and times. We present a simple analytic model for the runaway based upon 
a combination of stellar mixing-length theory and recent advances in understanding Rayleigh- 
Benard convection. The convective flow just prior to runaway is likely to have a strong dipolar 
component, though higher multipoles may contribute appreciably at the very high Rayleigh 
number (10^^) appropriate to the white dwarf core. A likely outcome is multi-point ignition with 
an exponentially increasing number of ignition points during the few tenths of a second that it 
takes the runaway to develop. The first sparks ignite approximately 150 - 200 km off center, 
followed by ignition at smaller radii. Rotation may be important to break the dipole asymmetry 
of the ignition and give a healthy explosion. 



Subject headings: Supernovae, hydrodynamics 

1. INTRODUCTION 

Despite forty years of study (Hoyle & Fowler 
1960), the mechanism whereby a degenerate 
carbon-oxygen white dwarf explodes, producing 
a Type la supernova (SN la), remains poorly un- 
derstood (for a recent review see Hillebrandt & 
Niemeyer 2000). Early calculations assumed that 
central carbon ignition would lead to a detonation 
(Arnett 1969) that would incinerate the star en- 
tirely to iron. This proved inconsistent both with 
observations of features in the supernova spec- 
trum from intermediate mass elements and with 



detailed calculations of isotopic nucleosynthesis. 
Nowadays it is understood that prompt detona- 
tion does not occur because the core at ignition 
time is insufficiently isothermal (Woosley 1992). 

Attention in recent years has thus focused on 
deflagrations, subsonic burning fronts in which 
pressure equilibrium is maintained across the 
burning interface. Though it is controversial 
whether the deflagration will later make a transi- 
tion to a detonation (Niemeyer & Woosley 1997; 
Khokhlov, Gran, & Wheeler 1997; Niemeyer 
1999), it is universally assumed that the runaway 
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begins as a deflagration (Nomoto, Sugimoto, & 
Neo 1976). 

There is less certainty, though, regarding ex- 
actly how and where the runaway is ignited. 
Most one-dimensional (ID) calculations, because 
of their imposed symmetry, obtain ignition at the 
center of the star. Many 2D studies have also as- 
sumed central ignition, largely as a matter of con- 
venience. However recent work (Niemeyer, Hille- 
brandt, & Woosley 1996; Rcincckc, Hillcbrandt, 
& Niemeyer 1999, 2002ab) has highhghted the 
sensitivity of the supernova outcome to precisely 
how the runaway is initiated - at the center or at 
one or more points off-center. These results may 
be summarized as showing that, in the absence of 
detonation, multi-point spherically-symmetric off- 
center ignition gives the most robust explosions, 
central single-point ignition gives weaker ones, and 
single-point off-center ignition gives such a weak 
explosion that it fails to unbind the white dwarf 
on the first attempt. 

It has been recognized for some time that the 
critical circiimstances affecting the ignition are de- 
termined when the convecting white dwarf core 
reaches a density of 2 — 3 x 10^ g cm~^ and a tem- 
perature of about 7 X 10® K (Nomoto et al 1984; 
Woosley & Weaver 1986). For these conditions, 
the time scale for the increase of nuclear energy 
generation becomes comparable to the convective 
turnover time, both of order 10 to 100 s. By the 
time any fiuid element reaches 10^ K, burnings has 
become quicker even than the time it takes a sound 
wave to cross a pressure scale height and, for all 
practical purposes, carbon burns instantly to iron. 
The surface of this fluid clement then becomes a 
"flame" , with a well determined speed (Timmes & 
Woosley 1992) and a buoyancy given by its density 
decrement, Ap/p ^ 15%. The explosion is born. 

Woosley (1990) first suggested that the much 
smaller buoyancy of convective fluid elements 
would still lead to appreciable radial motion, even 
as the convection decoupled from the burning. 
Thus ignition would occur off-center at one or 
more points already moving rapidly outwards. 
This speculation was rendered more quantitative 
by Garcia-Senz & Woosley (1995), who, using 
a simple parameterized description of "burning 
floating bubbles" , estimated a typical ignition ra- 
dius ~200 km. They also noted that off-center 
ignition would help alleviate a chronic overpro- 



duction of neutron-rich isotopes in Type la super- 
novae, since most of the burning would take place 
farther out at lower density than in centrally ig- 
nited models. 

Woosley (2001) estimated the convection speeds 
(~100 km s~^) and temperature fluctuations 
(AT/T ~0.3 - 3%) in the convective core at the 
time of runaway and, comparing them to defla- 
gration flame speeds and accelerations due to 
off-center burning, all of which are comparable, 
concluded that multi-point off-center ignition was 
probable. He also pointed out that the exact 
number and location of points was sensitive to 
small variations in the initial conditions, thus in- 
troducing some degree of chaos in the outcome. 
Since such observables as the kinetic energy, nickel 
mass and peak luminosity are sensitive to how the 
star ignites. Type la supernovae, starting from 
nearly identical initial conditions, will always ex- 
hibit some irreducible diversity. 

The conclusion that multi-point off-center ig- 
nition is likely was challenged by Hoflich & Stein 
(2002). Using a 2D implicit hydrodynamics code 
to follow the last few hours of the white dwarf 
runaway, they found no evidence for multiple spot 
or strong off-center ignition. Instead their model 
ignited at about 30 km, virtually at the center, 
and only once. Ignition was "induced by compres- 
sional heat" . As we shall see, their results, though 
a major computational advance, may have been 
influenced by attempting to model a 3D, spherical 
problem while carrying only a fraction of the solid 
angle on a 2D grid. The resolution and Reynolds 
number may also have been too low to see multi- 
point ignition in a simulation that was, at best, 
mildly turbulent and included only a small frac- 
tion of the fluctuation distribution function for the 
temperature. 

In this paper and its companion (Kuhlen, 
Woosley, & Glatzmaier 2003a; Paper II), we ex- 
plore the ignition of a nuclear runaways in Chan- 
drasekhar mass white dwarfs using two different 
approaches - an analytic model (this paper) and 
a 3D anclastic numerical model. In the analytic 
case, we are influenced by recent developments 
in understanding convection in Rayleigh-Benard 
experiments. These experiments show that the 
qualitative character of convection may change 
markedly depending upon the Rayleigh number. 
Even at high Rayleigh number, a persistent "roll" 
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dominates the flow pattern in Rayleigh-Benard 
convection. The analogue to this in a sphere is a 
dipole flow pattern. Recent numerical simulations 
(Paper II) suggest that such a large scale flow is 
also present in (non-rotating) convective stars. If 
so, it affects the mechanics of white dwarf ignition 
in a major way that can only be seen by carrying 
the entire sphere in the calculation. 

2. AN ANALYTIC MODEL OF THE 
RUNAWAY 

The flnal stages of the carbon runaway, wherein 
roughly 1.1 Mq of the core becomes convective, 
goes on for well over a century. By the time the 
central temperature reaches Tg — T/IO^K = 7 and 
pg = p/lO^g cm~^ = 2, the typical time scale for 
convection to go a pressure scale height, about 450 
km, is ~10 s, and has become comparable to the 
nuclear time scale. 

2.1. Nuclecir energy generation rate and 
time scale 

Nuclear energy generation during carbon igni- 
tion is given entirely by the (highly screened) fu- 
sion of two ^^C nuclei to form, chiefly, ^°Ne and 
^^Mg. The approximate energy generation rate 

(assuming carbon burns to a mixture of 3 parts 
2°Ne and one part ^^Mg; Woosley 1986) is 



6.7 X 10^^ X^(^^C) p9 F,^ Ai2,i2 



erg g ^ s \ 



(1) 



where A12.12 is the carbon fusion reaction rate 
(Caughlan & Fowler, 1988), Fsc is the electron 
screening function, and X(^^C) is the mass frac- 
tion of carbon. For a range of temperatures, Tg = 
6 - 8, 



A 



12,12 



7.6 X 10"^^ ( — 



30 



(2) 



The electron screening function (Alastuey & Jan- 
covici, 1978) is given by (pg = 1 - 3; Tg = 6 - 
8) 
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'P9 



2.3 



(3) 



so that the energy generation rate for a composi- 
tion of 50% carbon, 50% oxygen is 



5nuc « 2.8x10^3 



23 

(f) 



3.3 



erg g ^ s ^ 



(4) 



The specific energy available is 

qnnc = 4.0 x 10^^ XC^C) erg g-^. (5) 

The specific heat at constant pressure, which is 
required to estimate the nuclear time scale, is (e.g., 
Chiu 1968) 
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+ erg g-i (10«K)-\ 

P9 



(6) 



where A is the mean atomic mass number (13.7 
for 50% carbon, 50% oxygen by mass); m is the 
mass of the electron, x = pp/ mc is related to the 
mass density by 9.74 x IQ^peX^ = p; is the elec- 
tron mole munbcr, here 0.5; and the other symbols 
have their usual meanings. For the conditions of 
interest, e.g., Tg = 7 and pg = 2, the heat ca- 
pacity of the radiation field is negligible and the 
ions and electrons together provide 1.4 x 10^^ erg 
g^^ (10® K)~^. This relatively small heat capacity 
makes the carbon highly incendiary in the sense 
that a small amount of burning raises the temper- 
ature considerably. For lower densities, however, 
P9 ~ 0.01, the heat capacity of the radiation field 
becomes important, and this is what finally keeps 
the star from burning entirely to iron. 

For central temperatures near Tofi = 7, The 
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nuclear time under these conditions is 
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This is the time for the energy generation in an 

isolated region to increase from its starting value 
at temperature, Tg, to such high values that the 
reactions are virtually instantaneous. While con- 
vection remains efficient, this time is lengthened 
in the star by a factor of approximately 50 (§ 4.1). 

2.2. Luminosity 

The long convective episode in the pre-cxplosivc 
star establishes an overall adiabatic temperature 
gradient in the central regions. Using this condi- 
tion, the known density structure, energy gener- 
ation, Eq.(4), and assuming hydrostatic equilib- 
rium, one can estimate the luminosity. 

Because of the extreme sensitivity of the en- 
ergy generation to temperature, the luminosity 
will originate from a small fraction of the mass 
justifying a first order polytropic extrapolation (n 
= 3) of the central conditions. The equation of 
state is that of a relativistically degenerate gas. 



P = Kp^'^ 



= 1.24 X 10^^ (^yj dyne cm 



(8) 



and the polytropic radius parameter, 

( K \' 
a = I - 



= 385 km 



2/3 



Po,9 



1/3 



(9) 



with po^g the central density in 10^ g cm~^. Defin- 
ing a dimensionless radius, C = r/a, the mass in- 
terior to radius r is given by (^^ ^ 1), 

M{r) « f r-Vo(l - ^C'), (10) 



and the density at radius r is 

p{r) ^ Po{i - 



(11) 



Combining the equation of hydrostatic equilibrium 

dP _ GM{r)p{r) 

dr ~ r2 ' ^^^^ 

with the condition for an adiabatic temperature 
gradient 

f = (i-i/r.)(T/P)f, (13) 

one obtains, for r2 « 1.7, the variation of tem- 
perature with radius near the center of the white 
dwarf (Woosley 1990) 

r2-l 2TrGpy.h 

T2 5Po 

0.0185 (^)'^'Vir? 

'(14) 

where /i « (1 — j^C^) is a correction, near unity, 
for the density gradient. This equation also de- 
scribes the temperature evolution of any adiabat- 
ically expanding (or contracting) fluid element as 
it moves in a region near the star's center. 

One may then integrate Eq.(14) to obtain the 
luminosity as a function of central temperature, 



T(r) ^ Toil 



To 1 



L = 4tt j SnncPr^ 



dr 



7.0X10^4 g-l (^)4.3(Zk)23j^ 



7 



(15) 



where / is the integral, 

1 = jr^l - hr^.f^hdrr, (16) 

with b = 0.0185(p9/2)V3/i and /2 ~ (1 - K^)^'^- 
This integral can be evaluated numerically to give 
/ = 0.98 and 0.65 for po,g = 2 and 3 respectively. 
The result for the luminosity is valid to better than 
20%. 

One can also estimate the size of the energy gen- 
erating region by calculating the radius where L 
reaches one-half its value, 140 and 120 km respec- 
tively for pQ = 2 and 3. That is, approximately 
one-half of the luminosity of the star is generated 
in its inner 130 km. 
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2.3. Results from mixing-length theory 

As the runaway proceeds, the central tempera- 
ture rises and, along with it, the luminosity. En- 
ergy is transported by convection and dissipated 
by expansion, i.e., PdV work against gravity, and 
neutrinos, both of which occur chiefly outside the 
energy-generating central core. 

Conditions near the end of this ramp up can be 
estimated using mixing- length theory (e.g., Clay- 
ton 1983; Lantz & Fan 1999). In particular, the 
heat flux at radius r is given by 



= I ^;rms Cp (/AVT) 
« ^ Vrms CP (AT) 



(17) 



where L(r) is the luminosity at radius r, p is 
the density, Wrms, the average convective velocity 
there, cp, the heat capacity, Eq.(6), /, a character- 
istic size for the convection region, and AT, the 
temperature change across this region in excess of 
the adiabatic value. The factor "i" accounts for 
the fact that heat is carried by the outward mov- 
ing fluid elements, while lower entropy elements 
return the mass. Density is assumed nearly con- 
stant in all elements at a given radius. 

Fluid elements are buoyant because, at con- 
stant pressure, their excess temperature is accom- 
panied by a deficiency in density. The logarithmic 
derivative of density with respect to temperature, 
at constant pressure, is given by 



^P = - 

T 
P 
T 
P 



dlnp 
dlnT 

Ap 

AT 



zap 

T 



(18) 



w 1.9 X 10" 

P 

for a range of temperatures near Ts = 7 and 
P9 = 2. The derivatives in the above equation 
were evaluated numerically using the equation of 
state in the Kepler code (Weaver, Zimmerman, k, 
Woosley 1978). For example, for pg = 2 and at 
Tg = 7 and 8 respectively, 5-p = 6.6 x 10"^ and 



7.7 X 10""^ respectively. The small value of this di- 
mensionless constant reflects the extreme degen- 
eracy of the gas, i.e., that a large temperature 
change is required to give a pressure change com- 
parable to that resulting from a small change in 
density. 

A typical convective velocity is then given by 



2<?Ap 



1/2 



(19) 



with 5, the local acceleration due to gravity, i.e., 
GM{r)/r'^, and Ap, the density variation corre- 
sponding to AT. Combining Eqs. (17), (18), and 
(19), one has the mixing-length-theory estimate 
for the typical convection speed. 



AgrSp (f>{r) 
pcpT 



AGSpL 
3cpT 



1/3 



1/3 



(20) 



Here L is the luminosity in erg s ^. So long as L 
is nearly constant, the result is insensitive to I. 
For L45 = i(200 km)/10^5 erg s'^, 



40 km s"^ ( 
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Nl/3 rl/3 
) ^4 



^45 



(21) 



Since L is zero at the stellar center, Vj-^s formally 
goes to zero there, but its actual value depends 
upon the pattern of fluid flow. As is discussed in 
§ 3.2, the velocity just outside the energy generat- 
ing core, say at 100 - 200 km, may be a character- 
istic of the entire convection region. 

Most of the temperature dependence of Wrms 
is contained in the L/T term. This depends 
on the central conditions roughly as T^^p^-^, so 
the velocity scales with the central temperature 
roughly as Tjp^'^. From Eq.(15) the luminosity 
for To, 8 = 7, pg = 2 is about 10*^ erg s^^, so 



40kms-i (f)^-^^)^. 



(22) 



That is, one expects speeds about 2.5 times faster, 
or 100 km s^^, at Tq.s ~ 8 compared with To,8 = 7. 

The corresponding characteristic temperature 
excess across the convection zone is 



AT = 



p^c^gl 5p 



1/3 



(23) 



5 



or equivalently, 
AT 



2glSp 
= 0.005 (1) 



2.1 



14 



(24) 



This is the change in temperature excess across 
the region, beyond the adiabatic vahic. It is 
much smaller than the total change in tempera- 
ture across the convection zone. 

The typical density variation, Ap/p, is Sp times 
this, or about 3.6 x 10~^. Because of the uncertain 
choice of I and the use of a single energy character- 
ized by AT to carry the flow, these estimates are 
probably only accurate to a factor of two. There 
may also be some weak dependence of AT on the 
Rayleigh number not included in our simple anal- 
ysis. 

As we shall see later, it is a characteristic of 
the higher temperature fluctuations that ignite the 
explosion that they move with a higher speed than 
the average, Vrms- Also at To^g = 7.5 — 8.0 which 
may be a more typical choice for the mean central 



temperature at ignition, i\ 



100 km s-i. It is 



interesting, and seemingly an accident of nature, 
that, to about a factor of two, this is the same 
as the conductive laminar flame speed (Timmes & 
Woosley 1992), 



76 km s" 



-1 0-805 



2 7 



X(i^C) 
0.5 



(25) 

Partly because of this coincidence there is a depen- 
dence of the outcome on the carbon mass fraction 
in the white dwarf interior. 

2.4. Some characteristic measures 

The opacity, given by electron conduction 
(Timmes 2000 and 2002, private communication), 
is 

1.4 X 2.2 



5„™2 „-l / 2 \ ' /Ts 



Kcond » 2.7 X 10-^cm^ g-^ i^-j 

Using this to get the conductivity, 

4acT3 
o- = , 

PK 

for typical conditions, a w 3 x 10^^ erg cm~^ 
s~^. The viscosity, r], is given (Nandkumar & 



(26) 



(27) 



Pethick 1984) by 

1.9x109 /p9\ . 



■q 



(f) 



cm ^ s ^, 



(28) 



where I2 ~ 0.5 and Z « 7. Hence 77 w 10^ g cm"-*^ 
s~^ . From these one can calculate a Rayleigh num- 
ber (the so called "dimensionless temperature gra- 
dient" in convection; the ratio of buoyancy forces 
to diffusion forces) 



Ra 



g l-^ p^ CP 6p AT 



Trja 



(29) 



10 



,25 



where we have assumed that the appropriate AT 
is approximately (AVT) Z, Eq. (17), from mixing 
length theory. One can also estimate a Prandtl 
number (ratio of momentum transport to heat 
conduction) 



Pr = 



cprj 
a 

4 X 10"^, 



(30) 



a Reynolds number (ratio of inertial forces to vis- 
cous forces). 



Re = 



lol^ 



(31) 



and a Kolmogorov length (where turbulence dissi- 
pates). 



» 3 X 10-"* cm. 



(32) 



This combination of Rayleigh and Prandtl num- 
bers is well beyond the limits of what can be stud- 
ied on the Earth by either experiment or simula- 
tion. 

One can also estimate the Nusselt number (to- 
tal rate of heat transfer compared with conduc- 
tion). 



Nu 



>l 



a AT 

_ P Vrms I Cp 

2a 

w 3 X 10" 

which will be relevant in later discussions. 



(33) 
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3. LESSONS FROM RAYLEIGH-BENARD 
CONVECTION 

3.1. Recent Developments 

There exists a rich hterature of experiments and 
simulations that study the convection of matter 
between a hot and a cold plate. An important is- 
sue is how the relation between heat transport and 
the Raylcigh number, called Nu — Ra scaling, itself 
varies with Ra (e.g., Castaing et al. 1989; Gross- 
man & Lohse 2000; Kadanoff 2001). Another im- 
portant issue is how the form of the probability 
distribution finiction for temperature fluctuations 
(PDF) scales with Ra and Pr. Is it exponential 
or Gaussian? 

Below about Ra ~ 10^, it is thought that flow 

in a Raylcigh-Benard cell has not become com- 
pletely chaotic, but above 10* the regime of so 
called "hard turbulence" is encountered. Char- 
acteristics of hard turbulence include (Xia & Qiu 
1997): 1) A Nusselt number (Nu) that scales as 
Ra^/"^ (Castaing et al. 1989) instead of the classi- 
cally expected Ra^/^; 2) a PDF that IS exponen- 
tial, not Gaussian as might have been expected 
based on the central limit theorem (Kolmogorov 
1962); and 3) coherent large-scale circulation. Ex- 
periments by Niemela et al. (2000) confirm these 
properties of hard turbulence across a large range 
of Rayleigh numbers, 10^ - 10^^. Numerical calcu- 
lations by Rogers, Glatzmaier, & Woosley (2003) 
confirm that at least the Ra-Nu scaling character- 
istics of hard turbulence persist when the flow is 
compressible with density stratification. 

However, our white dwarf problem is character- 
ized by Ra still 10* larger than studied by Niemela 
et al. More importantly perhaps, we are interested 
in convection without hard boundaries. Energy is 
dissipated in the star by expansion and neutrinos, 
both volumetric losses, not conductivity to a cold 
boundary. Heating also occurs over an extended 
region making for a very thick "wall zone" and a 
range of temperature fiuctuations that is not lim- 
ited by the temperature of some hot plate. In the 
la, the flow can also pass right through the core of 
the burning region and out the other side, allow- 
ing more complete mixing between the core and 
the convective region than is possible in Rayleigh- 
Benard. 

A hint of what may lie ahead at very high val- 



ues of Rayleigh number, comes from studies by 
Kraichnan (1962), and is referred to in the lit- 
erature as the "ultimate" or "Kraichnan" regime 
of convection (Grossman & Lohse 2000; Kadanoff 
2001) . In situations where the effects of walls 
are suppressed, this transition may occur at much 
more modest values of i?a, even 10^ (Lohse & 
Toschi 2003). In this regime, it is thought that 
Nu oc (RaPr)^^'^ and that the large scale "rolls" 
seen in Rayleigh-Benard convection at lower Ra 
may give way to a more fragmented, chaotic pat- 
tern. It is unknown whether the distribution of 
temperature fluctuations in this regime is Gaus- 
sian or exponential. 

It is noteworthy, however, that the Nu — Ra 
scaling in this ultimate regime recovers the same 
simple scaling for AT found from mixing length, 
i.e., Eq.(23). That is, Eqs.(29), (30), and (33) plus 
the condition 

Nu ~ Rn}/^Pr^/^ (34) 

implies that 

/ y \ 1/3 

which to a factor of order unity is Eq.(23). This 
correspondence does not exist for any other scal- 
ings between Nu and i?a, which in general would 
leave some residual dependence of AT on the con- 
ductivity or viscosity. The mixing length approx- 
imation is, apparently, equivalent to Kraichnan 
scaling. Might the other properties of near in- 
finite Ra number circulation in Rayleigh-Benard 
cells also be relevant? 

3.2. Flow Patterns 

To a large degree, the velocities, temperature 
fluctuations, and, ultimately, the ignition process 
depend upon the circulating pattern assumed by 
the major flows within the core. We consider two 
representative cases - the "isotropic model" and 
the "dipole model" . The actual solution may have 
aspects of both. 

In the isotropic model, there is no preferred di- 
rection. Matter enters the burning region from all 
angles, is heated, reverses its direction and flows 
out. In the perfectly symmetric model, matter 
at the center is at rest, but this ideal state is 
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never realized. Small imbalances in flow will re- 
sult in overshooting, first in one direction, then 
another, so that the velocity vector at the origin 
varies in random way, being zero only on the av- 
erage. A slice at constant radius of the convective 
flow might resemble the solar photosphere, except 
that the large entropy contrast between upward 
moving plumes and downflows would be absent. 
For high Ra, the characteristic radial speed near 
the center would be less than, but perhaps not 
much less than that say a hundred km out. In this 
sense, the flow pattern might resemble the radial 
equivalent of Raylcigh-Benard in the "ultimate" or 
Kraichnan regime. Though not yet studied in the 
laboratory because of its inaccessibly high Ra, it 
is hypothesized that in this regime plumes would 
travel from one plate to another in near ballistic 
fashion with the large scale circulation suppressed 
(KadanofT 2001). 

One might expect that in the absence of rota- 
tion there can be no preferred orientation, that 
is, this isotropic model would be the only physi- 
cal one. However, numerical experiments (Kuhlen, 
Woosley, fc Glatzmaier 2003b; Paper 2) show that 
the convective flow in either a sphere or thick shell 
often takes on a dipole character. Matter flows in 
from one side, is heated in the central region, and 
flows out the other, like a jet engine. The center, 
far from being a point of stagnation, is character- 
ized by the same high velocities found farther out 
in the convection zone. This is the stellar analogue 
of the large scale circulation, or "rolls", seen in 
Rayleigh-Benard convection. Similar dipole flows 
have been seen in three-dimensional studies of 
red giant convection by Woodward and colleagues 
(Porter, Anderson, & Woodward 1997; Porter, 
Woodward, & Jacobs 2000; Woodward, Porter, & 
Jacobs 2002, 2003). 

This dipole circulation is an example of spon- 
taneous symmetry breaking. In the absence of 
rotation, the dipole picks an arbitrary axis, de- 
termined ultimately by tiny perturbations in the 
initial model, and maintains it for many convec- 
tive turnover times. Over very long periods, the 
orientation of the dipole may vary due to the occa- 
sional large "intermittent" occurrence. Also, the 
average over many calculations with random start- 
ing conditions would give no preferred angle, so, 
in a sense, the model is still isotropic. But dur- 
ing the time the runaway develops, the main flow 



may be highly directional and this has important 
implications for the explosion that follows. 

4. IGNITION 

The runaway first commences when the tem- 
perature of the hottest fluctuation, Tmax, leads to 
nuclear heating faster than the adiabatic cooling 
that occurs when a blob crosses the burning re- 
gion. Ignition will thus occur when the integral 
along a convective path. 



dr 



diverges. From Eq. (14) 



Sn 



dT 
dr 



Cp Ur, 



dr. 



exp 



-0.037Te(^|)2/3r7. 



(36) 



(37) 



It remains to specify the velocity and some dis- 
tribution of starting temperatures at the center of 
the star. As we shall see shortly (§ 4.3), the aver- 
age central temperature at ignition is in the range 
To,& = 7.5 - 8.0. From Eq. (19), v,^^ « 60 - 100 
km s~^. Taking 80 km s~^ as representative, we 
find that the first runaway will occur at a radius 
of over 300 km for a fluctuation hat started with 
a central temperature Tmax,8 = 8.623.... 

This estimate is a little large. The initial per- 
turbations in the star's center cannot have their 
temperature specified to arbitrary accuracy. As 
we shall sec in § 4.3, the typical variation in the 
temperature of the hottest few points in the core is 
a fraction of AT, Eq.(23), about 1%. Thus instead 
of Tmax.s = 8.623 one should realistically consider 
temperatures in the range Ti„ax,8 = 8.5 — 8.7. This 
implies typical ignition radii in the range 150 - 200 
km, though larger values are still possible in rare 
cases. 

The hottest fluctuations probably begin with 
(though do not necessarily end up with) radial ve- 
locities lower than the average. Longer residency 
in the burning region leads to higher temperature 
fluctuations. Doing the same calculation for an 
assumed average radial speed of 20 km s~^, one 
obtains ignition for hot fluctuations in the range 
Tlnax.s = 8.0 — 8.2 at radii again near 200 km. It is 
important for the issue of multi-point ignition that 
the igniting fluctuation takes from one to several 
seconds to reach its ignition radius (§ 4.4). 
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This maximum fluctuation temperature at run- 
away, Tmax, is to be distinguished from the mean 
central temperature, To, in that it reflects condi- 
tions in a small, atypical fluid clement out on the 
tail of the probability distribution function (PDF) . 
Just how far out, is a critical issue. 

4.1. The probability density function for 
fluctuations 

In order to address this question, one needs to 
assess the probability that a given high tempera- 
ture fluctuation will occur. Ignition does not hap- 
pen when the average core temperature reaches a 
flash point, but when the hottest persistent fluc- 
tuation does. 

We assume a simple model that should cap- 
ture the essence of the real runaway. Divide the 
star into two regions: the burning core with mass 
Mcore, and a larger reservoir, Mconv, of cooler mat- 
ter to which it is convectively coupled. Because of 
the finite speed of convection, Mconv shrinks dur- 
ing the final stages of the runaway while Mcore re- 
mains nearly constant, but we examine conditions 
during the last few seconds before the explosion 
when this ratio is approximately constant. 

Matter passing through Mcore gets a tempera- 
ture increment, AT, due to nuclear burning. The 
magnitude of this increment is sufficient, when 
mixed with Mconv, to heat the entire region by 
an amount AT during a convective mixing time. 
That is, each time the convection zone circulates 
once, one must generate an amount of energy equal 
to its adiabatic excess, AVT which is lost. Since 
Mconv 3> Mcore) the time spent by a fluid element 
in the energy generating core, is much less than 
the convective mixing time, by approximately the 

ratio Mcore /Mconv 

An operational value for this ratio comes from 

comparing the time it takes matter at the center 
of the star to appreciably increase its temperature 
in the presence or absence of convective cooling. 
From a one-dimensional calculation that employs 
a time-dependent mixing-length model for convec- 
tion (the Kepler code. Weaver, Zimmerman, & 
Woosley 1978), for a range of central tempera- 
tures 7.0 < To, 8 ^ 7.5 this ratio is approximately 
50. Very crudely this is also the the ratio of the 
mass within one pressure scale height, about 0.3 
M0, to that inside 130 km, the one- half energy 



generating region, about 0.01 M©. Henceforth we 
adopt Mconv /Mcore = 50. The time that a fluid el- 
ement resides in Mcore is approximately its radius 
divided by the convective speed, 100 km s"-'^, 
or 1 second. The total convective mixing time for 
Mconv is thus ~50 seconds. 

During each mixing event in Mcore, a certain 
fraction of the material, / < 1, is exchanged. 
Given the efficient nature of convection and the 
turbulent nature of the core, / may be close 
to unity. A representative value 0.9</<0.99 
will be assumed. This implies that each mixing 
leaves behind between 1% and 10% of the mass 
in Mcore- This residual matter increases its tem- 
perature beyond the average for the core, which 
in fact changes very little. After n mixing events, 
each approximately 1 s in duration, the fraction 
of material that has increased its temperature by 
nAT is (1 — /)". Such a power law naturally gives 
rise to an exponential PDF for temperature fluc- 
tuations in the core, such that the probability per 
unit volume (or mass) for finding a fluid element 
with temperature, T, when the average is Tq, is 

EPDF = --^^^^A exp Ml - f){T - n)/AT] , 

(38) 

Here, the PDF has been normalized on the tem- 
perature domain (To, oo). This function gives the 
probability, per unit temperature, that a measure- 
ment at a random time and place will give a tem- 
perature, T. The integral of this function from 
T > To to infinity is the fraction of the mass that 
has temperature T or greater. Because of the sim- 
ple toy model assumed here, the PDF is exponen- 
tial (hence, "E" in EPDF), but one cannot exclude 
that it might instead be Gaussian. If the fluctua- 
tions are statistically independent of one another, 
the central limit theorem implies 

(39) 

where F < 1 is an uncertain factor, not much 
less than one, to be determined by experiment or 
simulation. 

Unfortunately, the PDF by itself does not say 
how the temperature excess is distributed - in a 
few large blobs or very many small ones. For 
this, one needs additional information or assump- 
tions. In particular, how many thermally discrete 
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regions exist in the inner 100 km where the energy 
is mostly generated? 

4.2. Persistence of fluctuations 

A minimum size for the thcrmaUy discrete cell 
is capable of crossing the burning region while re- 
taining its excess heat for a time, r, without cool- 
ing by conduction is Acondi given by 



\ 2 \ 2 

'^cond _ ^cond P 

D a 
I 

^ 1 s, 



(40) 



or Acond ~ 10 cm. Based upon this scale alone, 

the number of thermally discrete regions in A^core 
would be very large, ~ 10^*, but turbulence 
changes this by making the size of the thermally 
discrete region much larger.. 

In a turbulent fluid, large eddies help to gen- 
erate fluctuations which the small ones disrupt. 
Consequently, there is some cut off to the size 
spectrum of persistent fluctuations that depends 
on the intensity of turbulence in the medium. A 
relevant length scale is the distance over which 
the average temperature varies by a characteristic 
fluctuation scale. For a nearly adiabatic core. 



Aturb 



AT 
VTad' 



(41) 



For AT given by Eq. (24) and VT^d by Eq. (14), 
Aturb is about a km. 

Eddies smaller than this will tend to smooth 
out fluctuations larger than AT, while larger ones 
will just move them around and create new ones. 
The number of such regions inside 100 km that 
might survive against turbulent mixing is about 
10^. This is a more realistic estimate than the 
larger number given by conduction alone, and will 
be used in subsequent discussions. Fortunately the 
answer will only depend on the logarithm of this 
number. 

Of course, the turbulence is not really homo- 
geneous or isotropic. Porter & Woodward (2000), 
in three-dimensional numerical studies of stellar 
convection up to Ra ~ 10^^, find that the up- 
welling plumes, the ones that would carry the ig- 
nition sparks here, are considerably more laminar 
than the highly turbulent downdrafts. It is diffi- 
cult to quantify this result for the present problem 



of white dwarf ignition, but regions of a km or so 
might survive turbulent dispersal for the several 
seconds it takes them to run away. 

4.3. Exponentiation of ignition points 

Once a single point has ignited, will more follow 
or will the runaway only ignite once? When the 
winner crosses the finish line, how close behind are 
the second, third, etc. runners? 

We previously derived (§ 4) Tmax.s ~ 8.0—8.7 as 
the central temperature of the unusually hot spark 
that finally ignites the runaway. We shall adopt 
7max,8 = 8.5 in what follows, but similar results 
are obtained for other values in this range. At 
the time of interest, there will be a large number 
of points in the core with temperature close to, 
but less than Ts = 8.5. In fact the average central 
temperature at that time will be Tg^g = 8.5 — nAT 
where 

" = ind^r <*^' 

For N = (100 km/A)'' ~ 10*" and assuming an 
exponential PDF with / = 0.9 - 0.99, n = 3 - 6, 
which should be typical. Even if '--^ 10^^ and 
/ = 0.9, n — 18, an extreme upper bound. 

From Eq.(24) one can evaluate AT/T to de- 
termine the range of average central temperature 
when the runaway finally commences to be 



8.5 
7.7- 



- (3 to 6) AT 
7.9. 



(43) 



The exponentiation time scale, that time during 
which the number of ignition points will rise one 

e-fold, is the time it takes for the temperature near 
Tg = 8.5 - dT to rise by ST, where from Eq. (38), 

5T = AT/|/n(l-/)|. (44) 

That is, for / = 0.9 and an exponential PDF, 

cp 



'exp 



Sn 



■ST 



0.1 s 



7.7 

2^0,8 



23 



2 \^-^ /^Ts 



Po,9 



0.02 



(45) 



A PDF that is Gaussian(Eq. 39), or a larger value 
of / would give a shorter time scale. 

For reasonable choices of / and N, the answer 
is considerably less than the time it takes a spark 



10 



to traverse the convection region, about 1 s, and 
comparable to the time it takes for the supernova 
itself to expand and shut off the runaway (§ 4.4). 
Hence it seems that multi-point ignition is favored, 
though, within current uncertainties, not guaran- 
teed. The number of ignitions per second will be 



reduces the temperature by 



h{t) = exp(f/rexp)- 



(46) 



As noted previously, this answer is relatively in- 
sensitive to the uncertain value of N (Eq. 42). 
Numerical calculations to confirm this result will 
need sufficient resolution to follow the evolution of 
temperature fluctuations far out on the tail of the 
PDF, but are badly needed. 

4.4. The end of ignition 

Nuclear burning in the core continues to raise 
the temperature for a time, but once ignition has 
occurred at any point, the spread of the flame 
rapidly reduces the binding energy of the white 
dwarf. Regions not yet encroached by flame arc 
cooled by this expansion, shutting off the ignition. 
Prom that point onwards, the calculation is purely 
one of flame propagation. 

When does this occur? We carried out a simu- 
lation of a 1.38 M0 white dwarf at the onset of car- 
bon runaway (To^g = 7.6; po.g = 2.5) using the Ke- 
pler stellar evolution code (Weaver, Zimmerman, 
& Woosley 1978). At that time, the energy input 
from nuclear burning, i.e., the convective luminos- 
ity had reached 8 x 10^^ erg s~^, in good agreement 
with Eq.(15). We then switched off nuclear energy 
generation and took a tiny time step, allowing the 
model to expand adiabatically with the existing 
velocity structure. Prom this, the rate of adia- 
batic reduction in central density appropriate to 
this luminosity was determined numerically. 



dlnp 
dt 



= -2.2 X 10"'' s 



4 „-l 



(47) 



From the equation of state, the adiabatic relations 
for dlnp and dlnT are 



dlnP _ p _ 4 
dlnp ^ 3 



dlnT 



(48) 



0.43. 



dlnP T2 
So, expansion in the absence of nuclear burning 



dlnT 
dt 



= -1.2 X 10"^ s-^ 



(49) 



Note that this relation between 6T and 6p is very 
different than Eq.(18) which is evaluated for non- 
adiabatic expansion at constant pressure. 

This rate of cooling by expansion should be pro- 
portional to the change in net binding energy, and 
thus to total rate of energy deposition in the star. 
Once the flame forms, this occurs, off center, at 
a rapidly increasing rate. Because all motions are 
initially very subsonic, the whole star responds, 
even at its center, to energy deposited anywhere. 

In the middle of the star, steady nuclear burn- 
ing was actually raising the temperature at a rate 
that can also be determined numerically, with the 
energy generation turned back on. 



dlnT 
dt 



= 3.1 X 10-2 s-^ 



(50) 



This implies that when the energy generation by 
the "flame sheet" becomes 250 times greater than 
L = Sx 10"^^ erg s"^, i.e., 2 x iC^ erg s"^, expan- 
sion will start to quench the ignition. 

A flame, once born, delivers energy at a rate 



e = qnncAup 



(51) 



where qnnc is the energy released by carbon and 
oxygen fusing to iron, about 7 x lO^'' erg g~^, A 
is its surface area, and u is the average speed nor- 
mal to the surface bounded by A. Niemeyer et 

al (1996), find, for single point off-center ignition, 
u ~ 200 km s~^ after about 0.1 s when the burned 
region is about 10^^ cm^ (their Fig. 4), so that 
e - 3 x 10"^ erg s'^. 

The ignition process thus ceases about 0.1 s af- 
ter the first spark ignites off center. Other sparks, 
however, will still be in transit since it takes 1 s 
to emerge from the core (§ 4). These hot, localized 
burning regions, already in the process of running 
away, will be more difficult to extinguish. We es- 
timate that the total ignition process goes on for 
several tenths of a second. 

4.5. Angular distribution of the sparks 

Even though the exponentiation time scale in 
Eq.(45) is short and the number of potential ig- 
nition points large, there is some ambiguity in 
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the counting of discrete points. If the dominant 
flow is, for example, dipole in nature (§ 3.2), ig- 
nition will occur preferentially on one side of the 
star. The concept of multi-point ignition becomes 
blurred if those points are all in close proximity. 

Indeed there may be great difficulty getting a 
viable supernova explosion if all the ignition occurs 
on one side (Niemeyer, Hillebrandt, & Woosley 
1996). Unless a transition to detonation occurs or 
pulsational oscillations, it will be difficult to ever 
burn the other side. The explosion will then be 
sub-energetic and produce too little ^^Ni. 

We speculate that the symmetric large scale 
flow, the dipole term in the case of spherical con- 
vection, is broken at the high Rayleigh number 
characteristic of the white dwarf as it runs away. 
In this sense the dipole flow in stellar convec- 
tion is like the large cell-sized circulation seen in 
Rayleigh-Benard experiments. The flow is dom- 
inated by the largest symmetric scales. High 
Rayleigh number as in Eq. (29) may lead to di- 
minished importance of both. It is also possible 
that rotation will diminish the effect of the dipole 
(Kuhlen, Woosley, & Glatzmaier 2002). 

If so, the fluctuations flowing out of the core 
may be not only numerous and small, but nearly 
isotropic. Experiments and numerical simulations 
to test this speculation are needed, but will be 
difficult. 

4.6. Further evolution of sparks 

Our calculations suggest multi-point ignition 

roughly 150 - 200 km off center. This is a small 
amount of mass, at most a few hundredths of a 
solar mass. Is the distinction between r = and 
150 km important? 

Yes, it is. Once a spark burns to iron, its density 
contrast, rather than being Ap/p ~ 10~* will be 
15%. The radial acceleration resulting from this 
buoyancy will be 

" ^ [t) (52) 
w 8 X 10^ r7 cm s"^. 

A blob starting off center accelerates rapidly and 
does most of its burning far out in the supernova. 
In a tenth of a second it is already moving sev- 
eral hundred km s~^. One started dead center 
only expands initially at the laminar speed (or at 



a typical convcctivc speed, which is comparable) 
and, in the same 0.1 s, goes only about 10 km. 

5. CONCLUSIONS 

Based upon a simple analytic description of the 
runaway, we expect that the first ignition will oc- 
cur roughly 150 — 200 km off center when the 
hottest point in the center of the star reaches 
8.0 - 8.7 X 10* K (§ 4). The average central tem- 
perature of the star when these fluctuations are 
produced is To,s « 7.7 - 7.9 (§ 4.3). The orienta- 
tion of this point with respect to the star's center 
may be given by the existence of a strong dipole 
asymmetry to the macroscopic convective flow just 
prior to the nmaway (§ 3.2). 

The existence of this dipole component (§ 3.2), 
seen previously in numerical models of red-giant 
convection, in main sequence massive stars, and 
further explored in Paper II, is important in deter- 
mining the outcome of the explosion. Even given 
multiple ignition points, if they all happen on the 
same side of the star in close proximity, and if 
there is no subsequent detonation, a weak explo- 
sion will ensue. It may be, however, that quasi- 
symmetric flow is restored at higher values of Ra 
than have been simulated thus far. By analogy 
to Rayleigh-Benard convection, such a restoration 
of symmetry might be expected in the Kraichnan 
regime (§ 3.2). Alternatively, rotation may be in- 
voked to alter the dipole flow (Paper II). 

The convective speed, which gives a measure of 
the turbulent energy input on large length scales, 
is 50 to 100 km s"'^ throughout most of the mass 
of the white dwarf when it explodes. This, not 
the laminar flame speed gives a minimum rate at 
which the flame spreads, even in the absence of 
strong Rayleigh- Taylor instability. 

To carry the stellar luminosity at this speed, 
a superadiabatic excess, AT, is required across 
the burning region. Mixing- length theory and 
Nusselt- Rayleigh scaling give a common prediction 
for the characteristic super-adiabatic excess, AT, 
in the convection zone, i.e., Ra = {NuPry/^, only 
in the Kraichnan regime. All other scalings give 
a value for AT that depends upon the conductiv- 
ity or viscosity or both. Assuming the validity of 
(NuPrY^^ scaling, and based upon a simple toy 
model (§ 4.1), we show that the distribution of 
temperature fluctuations in the burning core may 
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be exponential with a characteristic temperature 
scale <AT. The hottest fluctuations are not pro- 
duced by compression as in Hoflich & Stein (2002), 
but by extended residency and burning near the 
center of the star. We give an expression for the 
expected PDF, Eq.(38). 

Once the first point has ignited, the delay until 
the second, third, etc. ignites is given by the time 
it takes a fluid element on the tail of the PDF to 
raise its temperature by an amount AT/| In (1— /)| 
where / is the fraction of the core's mass that gets 
exchanged during a mixing time scale. For rea- 
sonable assumptions, this delay is Tujn ^ 0.1 s, 
Eq.(45). This is small compared with the time 
it takes the flrst igniting spark to transverse the 
burning region (~ 1 s; § 4) and comparable to the 
time for the expansion of the supernova to shut off 
the ignition process (§ 4.4). The number of igni- 
tion events increases as approximately exp(f/Tjg„). 
Thus multi-point ignition is possible. 

However, an important unresolved issue is the 
angular distribution and mean separation of igni- 
tion points. This will be influenced by the dipole 
flow, but also by the propensity of hot regions to 
cluster. 

Given that the number of ignitions and their 
location is exponentially sensitive, by Eq.(45), to 
small variations in the conditions in the core at 
the time of runaway, one expects that SN la mod- 
els, starting from nearly identical circumstances 
will display chaotic variation in their properties, 
including their light curves. Factors that might 
cause a systematic variation in this dispersion in- 
clude the carbon-mass fraction (Eqs.l and 25) and 
the central density at ignition. Calculations to 
quantify these dependences will be difflcult, but 
should some day be undertaken given the impor- 
tance of Type la supernovae as standard candles 
for cosmology. 

Our analytic study can also serve to motivate 
and guide future numerical work. First, because 
the issue of the dipole nature of the flow is impor- 
tant, all calculations should carry the entire 47r 
solid angle of a sphere. Calculations that study 
just a wedge and assume spherical symmetry will 
miss an important aspect of the problem. Because 
of concerns that the dipole might be artificially 
created by the tendency of turbulence to cascade 
to large scales in 2D simulations, some calcula- 
tions, at least, will need to be done in 3D. 



A major goal of numerical simulation should 
be to ascertain the PDF for the temperature fluc- 
tuations. Is it exponential, as assumed here, or 
Gaussian? Is the characteristic scale given by AT 
in Eq.(44)? What is the operational value of /? 
What is the angular distribution and size distribu- 
tion of the fluctuations? The Rayleigh number of 
such studies needs to be as high as possible, flrst 
because Ra in the star is naturally unapproach- 
ably high, and second, to sec the dependencies of 
all the above answers, including the dipole flow, 
on Ra. In Paper II we shall address some of these 
questions. 
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